function zdvd = dvd(nleja,zl,tsc)

%Calculate divided differences for zl Leja points.
% nleja number of Leja points.
% tsc scaled time.
% zl Leja points.
 % zdvd the diveded differences.
% zfl Function of the Liouvillian to be aproximated.
% Remember!!! zdvd(k) is the coeff. of (z-zl_1)...(z-zl_(k-1)).


% f[zl(1)]
      zdvd(1)=exp(zl(1)*tsc);
% f[zl(1),zl(2)]
      zdvd(2)=(exp(zl(2)*tsc)-zdvd(1))/(zl(2)-zl(1)); 
      
       for jj=3:nleja
          zp=1.;
          zs=zdvd(1);

          for ii=1:jj-2
             zp=zp*(zl(jj)-zl(ii));
             zs=zs+zdvd(ii+1)*zp;
          end 
           
          zp = zp*(zl(jj)-zl(jj-1)) ;
         
          zdvd(jj)=(exp(zl(jj)*tsc)-zs)/zp;

       end 
   
       end 

